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Abstract 

The response of a coupled array of nonlinear oscillators to parametric excitation is calculated 
in the weak nonlinear limit using secular perturbation theory. Exact results for small arrays of 
oscillators are used to guide the analysis of the numerical integration of the model equations of 
motion for large arrays. The results provide a qualitative explanation for a recent experiment [Buks 
and Roukes, to appear in J. MEMS (2002)] involving a parametrically-excited micromechanical 
resonator array. Future experiments are suggested that could provide quantitative tests of the 
theoretical predictions. 
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I. MOTIVATION: NON-LINEARITY OF MEMS AND NEMS RESONATORS 



Recent technological advances have enabled the fabrication of mechanical resonators down 
to micrometer and even nanometer scales, with frequencies almost in reach of the GHz 
range.§'i Even though the main thrust in this field of research comes from the need to pro- 
duce smaller, lighter, faster, and more efficient electromechanical systems, there is new basic 
physics to be learned along the wayJi One particularly interesting aspect of the physical be- 
havior of micro- and nanoelectromechanical systems (MEMS and NEMS) is their nonlinear 
mechanical response at relatively small deviations from equilibrium. This nonlinear be- 
havior has been observed experimentally,!'! and also exploited to achieve mechanical signal 
amplification and mechanical noise squeezingi'i0 in single resonators. In addition, MEMS 
and NEMS facilitates the fabrication of large arrays of resonators, for which the coherent 
response might be useful for signal enhancement and noise reduction. It is important to 
understand the nonlinear behavior of MEMS and NEMS resonators in order to improve 
their future designs. At the same time, one can take advantage of these systems for the 
experimental study of nonlinear dynamics. 

This paper is motivated by a recent experiment by Buks and Roukesi (henceforth BR) 
who fabricated an array of 67 fully-suspended doubly-clamped micromechanical resonating 
beams. Each beam was 270/zmx lfimx0.25fim in size, and the distance between neighboring 
beams was 4/im. The substrate beneath the array was completely etched away, forming 
a suspended diffraction grating with optical access from both sides. All even-numbered 
beams were electrically connected to one electrode and all odd-numbered beams to a second 
electrode. This allowed the application of electrostatic forces to induce coupling between 
the beams. The system was driven parametrically by introducing an AC component to the 
potential difference between the even-numbered and odd-numbered beams. The collective 
response of the array, as a function of the driving frequency and the DC component of the 
potential difference, was measured using optical diffraction. The response that BR inferred 
from their measurement was surprising in that (i) instead of showing a band consisting of a 
sequence of resonance peaks at the 67 normal frequencies of the array, the typical response 
as the frequency was swept up showed a small number of wide peaks where the response 
gradually increased and very abruptly decreased; and (ii) the array responded at frequencies 
beyond the expected top edge of the band. 
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We show below that both of these effects are a direct result of the fact that the restoring 
forces acting on the resonators as well as the damping that they undergo are both nonlinear. 
In Section II we describe the simplest equations of motion that are required to model the 
nonlinear resonator array. In Section [ITT] we solve the response of a single nonlinear res- 
onator to parametric excitation at twice its resonance frequency using secular perturbation 
theory (for comparison, we solve in Appendix [A] the response to parametric excitation at 
the resonance frequency). In Section [TV] we use the same method to calculate the response 
of the coupled resonator array and obtain exact results for a few (two or three) resonators. 
Understanding the analytical results of these two sections allows us to interpret the results 
of Section where we numerically integrate the equations of motion for an array of 67 
resonators. Our results agree qualitatively with the observations of BR, explaining the two 
points mentioned above, but we suggest that further experiments be performed in order to 
test our theoretical calculations in a more quantitative fashion. 

II. EQUATIONS OF MOTION 

We seek the simplest set of equations of motion (EOM) that capture the important 
physical aspects of the array of coupled micromechanical beams. We first note that the 
normal frequencies of an individual beam are sufficiently separated that the frequency bands, 
formed by the coupling of the beams in the array, are well separated by gaps in which the 
system cannot respond. We therefore assume that we can treat the lowest band separately 
from all the others, so that each individual beam is oscillating strictly in its fundamental 
mode of vibration. Each beam can therefore be described by a single degree of freedom x n , 
giving its displacement from equilibrium. We neglect any inhomogeneities in the fabrication 
of the beams and assume that all beams are identical. BR have actually examined each beam 
individually and report that their beams have a fairly uniform distribution of resonance 
frequencies, with an average of ujb = 179.3kHz, and a standard deviation of 0.53kHz. There 
is a much larger variation in the quality factors of the beams prior to the application of 
electrostatic interaction between them, but this variation disappears when a small potential 
difference is introduced between the beams.@ 

The coordinates x n are all assumed small so that only terms to lowest order in x n , nec- 
essary to capture the physical behavior of the system, will be kept in the EOM. Two types 



3 



of forces act on the beams. The elastic restoring force of each beam and the electrostatic 
forces between the beams. Experiments done by Buks and Roukesi on single beams of the 
type used in the array show that their response is like that of a Duffing oscillator — an oscil- 
lator whose restoring force contains a term proportional to the cube of the displacement and 
makes the oscillator stiffer than it would be within the harmonic approximation. Assuming 
a symmetric restoring force, and therefore no term proportional to an even power of x n , and 
neglecting higher than cubic-order nonlinear corrections, the elastic force acting on the n th 
individual beam is 

F eiltic = -mu 2 B x n - max 3 n , (1) 

where m is the effective mass of a beam oscillating in its fundamental mode, whose frequency 
is u B . 

Even though the electrostatic force between two parallel charged wires decays only as 
1/r, for simplicity we consider only the attractive interactions between nearest-neighbor 
beams. Within this approximation each term in the EOM depends either on the variables 
x n , describing the displacement of an individual beam from its equilibrium position, or on the 
difference variables x n+ \ — x n , describing the relative displacements of a pair of neighboring 
beams. To keep the equations as simple as possible we restrict each type of nonlinear term 
in the EOM to depend either on x n or on x n+ \ — x n , depending on whether it is mostly 
influenced by the elastic forces of the beams or the electrostatic interaction between them, 
respectively. 

The cubic term in the expansion of the nearest-neighbor electrostatic interaction tends 
to pull the beams away from equilibrium, acting against the cubic term in the expansion of 
the elastic force in Eq. ([]]). Because, as we shall confirm later, the response of the array is 
consistent with having a cubic term which stiffens the beams, the elastic contribution to the 
cubic term is stronger than the electrostatic one. We therefore ignore the cubic as well as 
all higher terms in the electrostatic interaction, which we write as 

F dectric = ~mA 2 [l + Hcosu p t](x n+1 - 2x n + X n _L). (2) 

Note that the linear electrostatic force constant |mA 2 , which is modulated with a relative 
amplitude H -C 1, representing the DC and the AC components of the applied voltage, 
is positive, acting to soften the elastic restoring force. The factor of 1/2 is used with the 
difference variable for convenience. 
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Parametric excitation, as it appears in the bare Mathieu equation for a single oscillator 
of frequency ujq without all the additional terms that we have here, is an instability of 
the system that occurs whenever the drive frequency is around one of the special values 
uj p = 2ujQ/n, where n is an integer that labels the so-called instability tongues of the system 
(named after the tongue-shaped instability curves in the frequency-amplitude plane). We 
choose the parametric driving frequency uj p to be around twice some value uq within the 
array's band of normal frequencies. We are therefore exciting the system in its first (n — 1) 
instability tongue. Thus, 



where e is a small parameter. In the BR experiment the system was actually excited in its 
second instability tongue, i.e. u p was chosen around some frequency in the band. It turns 
out that the response at the second tongue, apart from a few differences, is quite similar 
to that of the first tongue. We therefore prefer to carry out full calculations only for the 
first tongue which is somewhat easier to handle, and just for comparison, we calculate in 
Appendix |A| the response of a single nonlinear oscillator, excited at its second tongue. 

There is good reason to believe that most of the dissipation in the coupled system is a 
result of the electrostatic interaction which causes currents to flow through the beams. This 
assumption is based on the observation of Buks and Roukesi that the quality factors greatly 
diminish as the DC component of the electrostatic potential is increased. We therefore make 
the simplifying approximation that dissipation occurs predominantly as a result of currents, 
all other dissipation mechanisms being relatively negligible. This approximation relieves 
us of the worry about the variation in the quality factors of the individual beams before 
application of the electrostatic potential. The dissipative forces in the EOM are therefore 
written with respect to the difference variable, 



where we have included a nonlinear contribution to the dissipation, of the same order as 
the nonlinear elastic force (JTJ) . When putting all the pieces together, we (a) divide out the 
effective mass m of a beam; (b) scale time t — > t/ujs so that all frequencies (including A) 
are measured in units of u>b] and (c) scale length x — > x/y/a to get rid of the dependence 



(3) 




+ -mu B ar][(x n+1 - x n ) 2 (x n+1 - x n ) - (x n - x Tl _i) 2 (x n - ir n _i)], (4) 
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on a. The equation of motion for the n beam becomes 

x n + x n + xl + ^A 2 [l + H cos(2u; + eQ)t] (x n+1 - 2x n + ar n _i) 

— 2"^^ n+1 — ^ n 

— — /j[(i n +i — x n ) 2 (x n+ i — x n ) — (x n — x n _i) 2 (x n — i n _x)] = 0. (5) 

In the following sections we shall solve these equations using secular perturbation theory. 
The physical parameter that allows us to use this approach is the linear damping coefficient 
which is assumed to be small. We therefore express it as T = 67, taking e to be our 
small expansion parameter. The parametric instability of the system then occurs for small 
driving amplitude near resonance, and if, in addition, we consider the system near the 
onset of the instability, we can assume that the effects of nonlinearity are small as well. 
Thus, for small displacements x n all the nontrivial physical effects, namely, the parametric 
excitation, the cubic elastic restoring force, and both the linear and the amplitude-dependent 
dissipation, all enter the EOM as perturbative corrections to the simple linear equations. 
All these perturbative terms can be chosen to enter the EOM in the same order of the 
small parameter e by taking the leading order in x n to be of order e 1//2 , and expressing 
A 2 H = eh. This ensures, as we shall confirm later on, that all the terms will contribute to 
the lowest-order solution we are seeking. The final form of the EOM is therefore 

X n + %n + %n + g [A 2 + eh COS(2^ + efl)t] {x n+ i - 2x n + X n -i) 

1 . 

— ip[{ x n+i - x n f(x ri+1 - x n ) - (x n - x n _i) 2 (i n - i n _i)] = 0. (6) 

As for boundary conditions, we follow the experiment of BR who had two additional fixed 
beams, identical to all the rest, at both ends of the array. This means that we define two 
extra variables and set them to zero, xq = xjy+i = 0. 



III. RESPONSE OF A SINGLE PARAMETRIC ALLY-DRIVEN NONLINEAR 
OSCILLATOR 

We begin by calculating the response of a single nonlinear oscillator to parametric excita- 
tion. Previous calculations of this problem exist in the literature0 (and references therein), 
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nevertheless, we solve it here as a precursor to the many-oscillator case, treated in the next 
section. The equation of motion (|6[) for the single-oscillator case becomes 

x + \uj 2 — e/icos(2c<j + eVt)t\x 

+ ejx + x 3 + r]x 2 x = 0, (7) 

where we choose ojq to be lo = \/l — A 2 , the resonance frequency of the beam in the harmonic 
approximation. The parametric excitation is performed around twice the actual resonance 
frequency of the oscillator. (In Appendix [A] we treat the case where the excitation is per- 
formed around the resonance frequency of the resonator.) 

We calculate the correction to linear response by using secular perturbation theoryEH 
Recalling that the motion of the oscillator away from equilibrium is on the order of e 1 / 2 we 
try a solution of the form 

x(t) = e 1/2 (A(T)e luJt + c.c.) + e 3/2 xi(t) + . . . (8) 

where T = et is a slow time variable, allowing the complex amplitude A(T) to vary slowly 
in time. The variation of A(T) gives us the extra freedom to eliminate secular terms and 
ensure that the perturbative correction Xi(t), as well as all higher-order corrections to the 
linear response, do not diverge (as they do if one uses naive perturbation theory). Using the 
relation 

: dA dA , . j _ * 

we calculate the time derivatives of the trial solution (^|) 

x = e 1/2 ([itoA + eA'}e iU + c.c.) + e 3/2 Xl (t) + ... (10a) 
x = e 1/2 ([-oo 2 A + 2iuoeA' + e 2 A"}e luJt + c.c.) + e 3/2 Xl (t) + ... (10b) 

Substituting these expressions back into the equation of motion (^)), and picking out all 
terms of order e 3 / 2 , we get the following equation for the first perturbative correction 

x x + aAi = - (2iuA'e iuJt + c.c.) + hcos[(2co + dl)t] (Ae iuJt + c.c.) 

- 7 (iuAe^ + c.c.) - (Ae iujt + c.c.) 3 - r) (Ae iujt + c.c.) 2 (iuAe iuJt + c.c.) . (11) 

The collection of terms proportional to e luJt on the right-hand side of Eq. ([TT]) , called 
the secular terms, act like a force, driving the simple harmonic oscillator on the left-hand 



side at its resonance frequency. The sum of all the secular terms must vanish so that the 
perturbative correction x\ (t) will not diverge. This gives us an equation for determining the 
slowly varying amplitude A(T). After expressing the cosine as a sum of exponentials we get 

d A h 

2iuy— - -^-A*e inT + iujA + 3\A\ 2 A + iur]\A\ 2 A = 0. (12) 

(J/J. _ 

Ignoring initial transients, and assuming that the nonlinear terms in the equation are suffi- 
cient to saturate the growth of the instability, we try a steady-state solution of the form 

A(T) = aS T . (13) 

The solution to the equation of motion (|7|) is therefore 

x(t) = e i/2( ae *(«+«n/2)* + c . c .) + o(e 3 / 2 ), (14) 

where we are not interested in the correction X\(t) of order e 3//2 , but rather in the fixed 
amplitude a of the lowest order term. This amplitude a can be any solution of the equation 

[(3|a| 2 - ojVL) + iufa + r]\a\ 2 )] a = ^a*, (15) 

obtained by substituting the steady-state solution (|13|) into the equation ( |T2"D of the secular 
terms. We immediately see that having no response (a = 0) is always a possible solution 
regardless of the excitation frequency fl We divide both sides of the last equation by juj 
and define the rescaled variables a = a/^fyoj, Vt = fj = cur], and h = h/2 r yu, in terms 

of which the equation for the fixed complex amplitude a becomes 

[(3|a| 2 i(l + f]\a\ 2 )] a = ha*. (16) 

Expressing a = \a\e % ^ we obtain, after taking the magnitude squared of both sides, the 
intensity \a\ 2 of the non-trivial response as all positive roots of the equation 

(3|a| 2 -n) 2 +(l + r / |a| 2 ) 2 = /i 2 . (17) 

This has the form of a distorted ellipse in the (Q, \a\ 2 ) plane, and a parabola in the (|a| 2 , h) 
plane. In addition, we obtain for the relative phase of the response 

i a* 1 1 + fj\a\ 2 

<p = - In — = — arctan , ( 18) 

v 2 a 2 3 a 2 
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FIG. 1: Response intensity \a\ 2 as a function of the frequency ti, for fixed amplitude h = 1.5. Solid 
curves are stable solutions; dashed curves are unstable solutions. Thin curves show the response 
without non- linear damping (fj = 0). Thick curves show the response for finite nonlinear damping 
(fj = 1). Dotted lines indicate the maximal response intensity |a|^ ax and the saddle-node frequency 



tation as a function of frequency f2, for fixed amplitude h = 1.5, in terms of the rescaled 
variables. Solid curves indicate stable solutions, and dashed curves are solutions that are 
unstable to small perturbations. Thin curves show the response without nonlinear damping 
(fj = 0) which grows indefinitely with frequency Cl and is therefore incompatible with the 
experimental observations of BR and the assumptions of our calculation. Thick curves show 
the response with finite nonlinear damping (77 = 1) . With finite fj there is a maximum value 
for the response |a|^ aa . = (h — l)/fj, and a maximum frequency, 



at which the stability of the solution changes (known as a saddle- node bifurcation). For 
frequencies above CIsn the only solution is the trivial one a = 0. These values are indicated 
by horizontal and vertical dotted lines in Figure [l]. 

The threshold for the instability of the trivial solution is easily calculated by setting 
a = in the expression ( |17D for the nontrivial solution. As seen in Figure [I], f° r a given 



In Figure [I] we plot the response intensity \a 



2 of a single oscillator to parametric exci- 
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FIG. 2: Threshold for instability plotted in the ($7, h) plane. The lower, long-dashed curve shows 
the threshold without any linear damping (7 = 0), which is zero on resonance. The upper curve 
shows the threshold with linear damping (7 ^ 0). The parameters for the upper curve are u = 1/2 
and 7 = 1 so that h = h. The threshold on resonance (0 = f2 = 0) is therefore h = h = 1. The 
solid and short-dashed regions of the upper curve indicate the so-called subcritical and supercritical 
branches of the instability, respectively. On the subcritical branch (J7 > fj/3) there will be hysteresis 
as h is varied, and on the supercritical branch (O < 77/3) there will not be any hysteresis. 

h the threshold is situated at €l = ± — 1. The threshold is plotted in Figure 0, in 
the (fl, h) plane. Note that the minimal amplitude needed for instability is obtained on 
resonance (O = 0) and its value is h — 1, or h — 2 / yu), so that it scales as the linear damping 
coefficient 7. 

Finally, in Figure |3] we plot the response intensity \a\ 2 of the oscillator as a function 
of amplitude h, for fixed frequency fl and finite nonlinear damping fj = 1. Again, solid 
curves indicate stable solutions, and dashed curves unstable solutions. Thick curves show 
the response for Q = 1, and thin curves show the response for Q = fj/3 and Cl = — 1. The 
intersection of the trivial and the nontrivial solutions, which corresponds to the instability 
threshold, occurs at h — vtW+l. For Q < fj/3 the nontrivial solution for \d\ 2 grows 
continuously for h above threshold and is stable. This is a supercritical bifurcation. On the 
other hand, for Q > fj/3 the bifurcation is subcritical — the nontrivial solution grows for h 
below threshold. This solution is unstable until the curve of \a\ 2 as a function of h bends 
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h SN [Eq. (20)] 
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FIG. 3: Response intensity |a| 2 as a function of the parametric modulation amplitude h for fixed 
frequency and finite nonlinear damping {fj = 1). Thick curves show the stable (solid curves) and 
unstable (dashed curves) response for = 1. Thin curves show the stable solutions for Cl = fj/3 
and Q = — 1, and demonstrate that hysteresis as h is varied is expected only for O, > fj/3. 

around at a saddle-node bifurcation at 



where the solution becomes stable and \d\ 2 is once more an increasing function of h. For 
amplitudes h <Jisn the only solution is the trivial one a = 0. 

Like the response of a forced Duffing oscillator, the response of a parametrically excited 
Duffing oscillator also exhibits hysteresis in a frequency scan. If the frequency Q starts out 
at negative values and is increased gradually with a fixed amplitude h, the response will 
gradually increase along the thick solid curve in Figure 0, until Cl reaches flsN and the 
response drops abruptly to zero. If the frequency is then decreased gradually, the response 
will remain zero until Q reaches the upper instability threshold \/h 2 — 1, and the response 
will jump abruptly to the thick solid curve above, and then gradually decrease to zero along 
this curve. A similar hysteretic behavior will be observed if the amplitude h is varied with 
a fixed frequency Cl > fj/3, as can be inferred from Figure |3]. 



1 + p 



(20) 
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IV. RESPONSE OF A PARAMETRIC ALLY-DRIVEN ARRAY OF NONLINEAR 
COUPLED OSCILLATORS SECULAR PERTURBATION THEORY 



Consider now the coupled array of nonlinear oscillators as described by the general 
EOM (|5|). We calculate its response to parametric excitation, again using secular per- 
turbation theory. We expand x n (t) as a sum of standing wave modes with slowly varying 
amplitudes 

N 

t (t) = e V2 (A m (T) sm(nq m )e^ 1 + c.c.) + e^ 2 x^\t) + ..., n = l...N. (21) 



m=l 



Recall that the boundary conditions are such that there are two additional fixed beams, 
labeled and N+l, exerting electrostatic forces on the first and the last beams of the array. 
With these boundary conditions (xo = xn+i = 0), the possible wave vectors q m are given by 

9w = iv+T' m = 1 --- N - ( 22 ) 



We substitute the trial solution (f2~TD into the EOM term by term. Up to order e 3 / 2 we 
have: 

x n = e 1/2 ^ sin(ng m ) ([-u 2 A m + 2iueA' m ]e iu)mt + c.c.) + e 3/2 x^(t), (23a) 



x n+1 - 2x n + x n _! = -4e 1/2 J2 sin2 ( y) sin(ng m ) (A m e iuJmt + c.c.) 

m 

+ e^(x[ 1 l 1 -2x^+x^ 1 ), (23b) 



^e 7 (x n+1 - 2x n + x n _i) = -2e 3/2 7 ^ sin 2 ) sm(nq m ) (iA m e luJmt + c.c.) , (23c) 



in 



x 3 = e 3 ^ 2 sin(ngj) sin(ng fc ) sin(ng/) 

x (Aje^ + c.c.) (A k e iuJkt + c.c.) (A t e^ lt + c.c.) 
e 3 / 2 

= — {sm[n(-g i + q k + %)] + sn\[n{qj - q k + ©)] 

+ sin[ra(g j + g fc - q t )) - sin[n(g i + q k + qi))} 
x { AjAkAe^^^ + SAjAkA^e 1 ^^-^ + c.c. } , (23d) 
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and 



yT] \{ x n+l ~ Xn) 2 (in+1 ~ in) 



{x n X n _i) {x n in— l)] 

Qj ■ Ik . 
— sm — sn_ 
2 2 2 



o 3/2 • Qj ■ Qk . Ql 



x ^sm 
+ sin 
+ sin 
+ sin 



j,k,l 

-qj + qk + qi 
2 

\ Qj~Qk + qf 
. 2 

"r gj + gfc-gt " 
. 2 
r<?j + % + Ql 



sm[n(-qj + q k + g;)] 
sin[n(g i - g fc + <#)] 
sin[n(g i + q k - q t )] 
sm[n(qj + q k + qi) 



(Aje^ + c.c.) (A k e iuJkt + c.c.) (iuiA^ + c.c.) .(23e) 



At order e 1 / 2 we simply get the linear dispersion relation, given by 



L), 



2 -l-2A 2 sin 2 (^), m = l...N. 



(24) 



At order e 3//2 we get N equations of the form 

(1) + + - A 2 (x^ii - 2x$ + x^l^j = ^ { mth secular term) e iuJmt + other terms, 



x. 



(25) 

where the left-hand sides are, again, coupled linear harmonic oscillators, with a dispersion 
relation given by fl2~3|). On the right-hand sides we have N secular terms which act to drive 
the coupled oscillators il 1 ' at their resonance frequencies. As we did for a single oscillator in 
section |IT|, here too we require that all the secular terms vanish so that the Xn remain finite, 
and thus obtain equations for the slowly varying amplitudes A m (T). To extract the equation 
for the m th amplitude A m (T) we make use of the orthogonality of the modes, multiplying all 
the terms by sin(ng m ) and summing over n. We also express all normal frequencies relative 
to the same reference frequency u , used to define the excitation frequency u p in Eq. (|3|), so 
that 



Mr, 



loq + eQ m . 



(26) 
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We find that the coefficient of the m th secular term, which is required to vanish, is given by 
- 2iuJ^ - 2^ m sin 2 + /^sin 2 (Sg) e^~^ T 



A A, A* J(^j+n k -Ui-U m )T A (l) 



j,k,l 

- 2 V { [2iu Jl A*A k A l e i ^ +nk+n ^ n ^ T - iu l A j A k A^ a ^"- ai -^ T ] 

3,k,l 

w A (2) . Qj . Qk . qi . q m \ n 

x A )ki; m sm 77 sm — sm-sm y| = 0, (27) 
where we have introduced two A functions, defined in terms of Kronecker deltas as: 



— $-j+k+l,m ~ 8-j+k+l,-m ~ 3- 


-j+k+l,2(N+l)-m 




"j—k+l,m Oj—k+l,—m ^j—k- 


H,2(JV+l)-m 




+ 3j+k-l,m ~ Oj+k-l,—m ~ fij+k- 


-l,2(N+l)-m 




~~ Oj+k+l,m + Sj+k+l,2(N+l)-m ~ 


' ^j+k+l,2{N+l)+rm 


(28a) 



and 



A (2) 
jkl;m 





0—j+k+l,m + 3-j+k+l,-m ~ 


-j+A;+i,2(JV+l)-m 




+ 


3j-k+l,m + 3j-k+l,-m ~~ Oj'-jfc- 


H,2(N+l)-m 




+ 


Oj+k—l,m $j+k— I,— m "j'+fc- 


-l,2(N+l)-m 




+ 


"j+fc+«,m — <5j+fc+«,2(7V+l)-m ~ 


' °j+k+l,2(N+l)+m- 


(28b) 



These A functions ensure the conservation of lattice momentum — the conservation of mo- 
mentum to within the non-uniqueness of the specification of the normal modes due to the 
fact that sin(ng m ) = sin(nq2k(N+i)±m) for any integer k. The first Kronecker delta in each 
line is a condition of direct momentum conservation, and the other two are the so-called 
umklapp conditions where only lattice momentum is conserved. 

As for the single oscillator, we again try a steady-state solution, this time of the form 

A m (T) = a m e<^) T , (29) 

so that the solutions to the EOM, after substitution of (|29"D into (]21|) , become 

,(*) = e 1/2 E ( a - sin(ng m )e^ + f> + c.c.) + 0(e 3 / 2 ), (30) 

m 

14 



X r 



where all modes are oscillating at half the parametric excitation frequency u p /2. 

As before, we are not interested in the corrections of order e 3 / 2 but only in the values of 
the fixed amplitudes a m as functions of all the parameters of the original EOM. Substituting 
the steady state solution Q2§D into the equations Q2"7D for the time- varying amplitudes A m (T), 
we obtain the required equations for the fixed complex amplitudes a m 



(Q - 2Q m )u m a m - 2i'yu m a m sin 2 + ha* m sin 2 - - ^ a j a kai^% m 



-2it] sin — 2^ oji [20^-afcaj - aja^ \ sin -± sm — sm — A^. m = 0. (31) 

i.fc,2 

We can change to rescaled variables as we did in the case of a single oscillator by dividing 
the equations for the amplitudes (|31|) by ('jujo) 3 / 2 and defining as before dj = a^j J^Uq, 
D, — fl/7, 77 = u;o?7, and fa = h/2juo, and in addition r m = cu m /u;o and 5 m = 2Q m /j. After 
doing so we obtain the rescaled equations 

(0 - 5 m )r m a m - 2ir m sin 2 a m + 2/isin 2 - - ^ n fika* ^%-,m 

j,k,l 

-2ifj sin y^n [2a*a fe a; - a y a fe a;] sin -j- sin y sin | A^. m = 0. (32) 

This is the main result of the perturbative calculation. We have managed to replace N 
coupled differential equations @ for the oscillator coordinates x n (t) by N coupled algebraic 
equations ( |3~TD for the time-independent mode amplitudes a m . All that remains, in order 
to obtain the overall collective response of the array as a function of the parameters of the 
original EOM, is to solve these coupled algebraic equations. 

Before doing so we should note the following general statements. First, one can easily 
verify that for a single oscillator (N = j = k = l = m = l), the general equation (|31|) reduces 
to the single-oscillator equation ( jl5]) , we derived in section [[II]. Next, one can also see that 
the trivial solution, a m = for all m, always satisfies the equations, though, as we have 
seen in the case of a single oscillator, it is not always a stable solution. Finally, one can also 
verify that whenever for a given m, A^ m .^ = A^ ■ = for all j 7^ m, then a single-mode 
solution exists with a m 7^ and cij = for all j 7^ m. These single-mode solutions have the 
elliptical shape of the single-oscillator solution given in Eq. (|17|), and satisfy the equation 

(-AM \n I 2 - q\ + (l + sin 2 *^A (2) n\a \ 2 Y - h 2 (W) 



4sin 4 (g m /2) V4 
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where for each solution we have set ujq = u m , so that 5 m = and r m = 1. Note that 
generically AmLm ; m = ^mmm;m = 3, except when umklapp conditions are satisfied. 

Additional solutions, involving more than a single mode, exist in general but are hard 
to obtain analytically. We calculate these multi-mode solutions below for the case of two 
and three oscillators by finding the roots of the coupled algebraic equations numerically. In 
Appendix |FJ we show the explicit sets of coupled mode-amplitude equations for these cases. 

In Fig. |] we show the solutions for the response intensity of two oscillators as a function 
of frequency, for a particular choice of the equation parameters. The top graph shows the 
square of the amplitude of the antisymmetric mode a,2, whereas the bottom graph shows the 
square of the amplitude of the symmetric mode d\. Solid curves indicate stable solutions and 
dashed curves indicate unstable solutions. The two elliptical single-mode solution branches, 
mentioned in the previous paragraph are easily spotted. These branches are labeled by Si 
and 5*2 (In Appendix [B], Equations ([B4D, we give the analytical expressions for these two 
solution branches). In addition, we find two double-mode solution branches, labeled D\ and 
D2, involving the excitation of both modes simultaneously. Note that the two branches of 
double-mode solutions intersect at a point where they switch their stability. 

With two oscillators we obtain regions in frequency where three stable solutions can exist. 
If all the stable solution branches are accessible experimentally then the observed effects of 
hysteresis might be more complex than in the simple case of a single oscillator. This is 
demonstrated in Fig. |5] where we compare our analytical solutions with a numerical integra- 
tion of the differential equations of motion (H) for two oscillators. The response intensity, 
plotted here, is the time and space averages of the square of the oscillator displacements 

< 34 > 

n=l 

where the angular brackets denote time average, and here N = 2. A solid curve shows 
the response intensity for frequency swept upwards, and a dashed curve shows the response 
intensity for frequency swept downwards. Small circles show the analytical response inten- 
sity, using the fact that / = 3( | o-i | 2 + |a 2 | 2 )/2, for the stable regions of the four solution 
branches shown in Fig. |]. With the analytical solution in the background, one can easily 
understand all the discontinuous jumps, as well as the hysteresis effects, that are obtained in 
the numerical solution of the equations of motion. Note the the Si branch is missed in the 
upwards frequency sweep and is only accessed by the system in the downwards frequency 
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FIG. 4: Two oscillators: Response intensity of as a function of frequency Cl, for a particular choice 
of the equation parameters. The top graph shows | c^2 1 2 , the bottom graph shows |ai| 2 . Solid curves 
indicate stable solutions and dashed curves indicate unstable solutions. The two elliptical single- 
mode solution branches [Equations ( p4a| ) and ( B4b| )] are labeled S\ and 52- The two double- mode 
solution branches are labeled D\ and Di- 
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FIG. 5: Hysteresis with two oscillators: Comparison of stable solutions, obtained analytically 
(small circles), with a numerical integration of the equations of motion (solid curve - frequency 
swept up; dashed curve - frequency swept down). Plotted is the averaged response intensity, defined 
in Eq. (|34|). Branch labels correspond to those in Fig. ||. 

sweep. One could trace the whole stable region of the Si branch by changing the sweep 
direction after jumping onto the branch, thereby climbing all the way up to the end of the 
Si branch and then falling onto the tip of the Di branch or to zero. 

In Fig. |6| we show the solutions for the response intensity of three oscillators as a function 
of frequency, for a particular choice of the equation parameters. The graphs show the squares 
of the amplitudes of the three different modes. Solid curves indicate stable solutions and 
dashed curves indicate unstable ones. For three oscillators there is only one elliptical single- 
mode solution branch, of the form of Eq. (^), whose exact analytical expression is given in 
Eq. (P8|). This branch is labeled by 5*2. In addition, we find a host of nontrivial multi-mode 
solution branches, including one that is disconnected from all other branches. We show 
these plots, not only to demonstrate that it is possible to obtain such solutions exactly, but 
also to emphasize the large number and nontrivial structure of the solution branches one 
finds, even for such a small number of oscillators. This can only serve as a hint for the 
multi-mode solutions one can expect to find when the number of oscillators is large, as in 
the BR experiment. 
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FIG. 6: Three oscillators: Response intensity of three oscillators as a function of frequency Cl, 
for a particular choice of the equation parameters. The graphs show the squares of the amplitudes 
of the three different modes. Solid curves indicate stable solutions and dashed curves indicate 
unstable ones. The only elliptical single-mode solution branch [Eq. (|B§|)] is labeled by S^- 
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FIG. 7: Response intensity as a function of the driving frequency to p (measured in units of the top 
band edge frequency Uq) for TV = 67 parametrically-driven oscillators (solid curve - frequency swept 
up; dashed curve - frequency swept down). The response intensity is defined as (x^) [Eq. fl34|)], 
with the bar denoting the average over the space index n, and the brackets the average over time. 
The parameters used are A 2 = 0.02, h = 0.16,7 = 0.004, rj = 6. 

V. RESPONSE OF PARAMETRICALLY-DRIVEN NONLINEAR COUPLED 
OSCILLATORS NUMERICAL INTEGRATION OF THE EQUATIONS 

The equations of motion @ were integrated numerically for an array of 67 oscillators, as in 
the BR experiment. The results for the response intensity ( |34"D as a function of parametric 
drive frequency u p (measured in units of the top band edge frequency ujq) are shown in 
Figure [7]. These results must be considered illustrative only, since many of the parameters 
of the experimental system are not known. The parameters used to construct the figure, 
A 2 = 0.02, h = 0.016, 7 = 0.004, rj = 6.0, were chosen using the insights gained from the two 
and three beam results. We should emphasize that the structure of the response branches 
depends very strongly on the equation parameters. 

A number of the important features of the experimental data are reproduced by these 
calculations. We concentrate on the solid curve in the figure, which is for frequency swept 
upwards, since this is the protocol that was used in the experiment. In particular, the 
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response intensity shows features that span a range of frequencies that is large compared 
with the mode spacing (which is about 0.0006 for the parameter used). The lowest frequency 
feature, from about u p /u! = 1.94 to uj p /u = 1.97, can be identified as the response to the 
parametric drive of a single mode at or very near the band edge at uj/ujq = 0.98, analogous 
to the one mode response shown in Figure |I| Furthermore, the variation of the response 
with frequency shows abrupt jumps, particularly on the high frequency side of the features 
as the frequency is raised. Finally, the response extends to frequencies higher than the band 
edge for the linear modes, which would give a response only up to uj p /u = 2.0. All of these 
features are understood now that we have seen the analytical solutions for small numbers of 
oscillators. In particular, the wide features compared with the mode spacing are explained 
by the simple fact the as the frequency is swept upwards a particular solution branch is 
followed as long as it remains stable. In the meantime many other stable solutions, that 
may be as close to each other as the mode spacing, are simply skipped over. 

Comparing the two traces in Figure [7] shows that the response for a downward frequency 
sweep is significantly different, with a less dramatic variation of the response. In particular, 
note that the downwards sweep was able to access additional stable solution branches that 
were missed in the upwards sweep. There is also no response above uj p /ujq = 2.0 in this case. 
This is because the zero displacement state is stable for uj p /ojq > 2.0, and the system will 
remain in this state as the frequency is lowered, unless a large enough disturbance kicks it 
onto another of the solution branches. The hysteresis on reversing the frequency sweep was 
not looked at in the first experiments, and it would be interesting to test this prediction in 
further experiments. 



VI. CONCLUSIONS 



We have calculated the response of nonlinear coupled oscillators to parametric excita- 
tion. Our calculations agree qualitatively with the experimental measurements of Buks and 
RoukeJl and explain the main features observed in the experiment. The abrupt drops in 
the response as the frequency is swept upwards, the response continuing beyond the upper 
edge of the frequency band, and the large size of the response features compared with the 
mode spacing, are all qualitatively explained. 

Nevertheless, we propose that a more systematic study be conducted on systems of cou- 
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pled nonlinear resonators so that our theoretical predictions could be tested more quanti- 
tatively. For example, successive measurements on systems containing only one, two, and 
three coupled resonators, that are made as identically as possible, could be used to extract 
the nonlinear parameters of the resonators. These could then be used to predict and explain 
the response of a large resonator array more quantitatively. 

Furthermore, we have demonstrated that as the number of oscillators is increased the 
number of the solution branches for the response of the system increases and the effects 
of hysteresis become more and more complicated. This suggests that the appropriate ex- 
perimental protocol for studying a system with many oscillators should be — in addition to 
the standard up-sweep and down-sweep in frequency — to change the direction of the sweep 
after every abrupt change in the response intensity. Such a protocol may provide more 
information about the response curve by accessing additional branches of the solution and 
fully tracing them out. 



APPENDIX A: PARAMETRIC EXCITATION OF A SINGLE OSCILLATOR AT 
ITS SECOND INSTABILITY TONGUE 



For a single nonlinear oscillator, like the one studied in section |TTT|, which is parametrically- 
excited at its second instability tongue, Eq. @ becomes 

x + [uo 2 - A 2 H cos(lu + efi)i] x 

+ e^x + x 3 + r\x 2 x = 0, (Al) 



where again lu = y/1 - A 2 is the resonance frequency of the oscillator in the harmonic 
approximation, but the parametric excitation is performed around uo and not around 2uo. 
In this case the scaling of A 2 H with respect to e needs to be redetermined. The technical 
reason for this is that if we naively take A 2 H = eh, as before, then the parametric driving 
term does not contribute to the order e 3//2 secular term which we use to find the response, 
and the order e 1//2 term in x becomes identically zero. 

The remedy for this situation is to let A 2 H scale like e p , with p < 1, so that there will 
be a non-secular correction to a; at a lower order than e 3 / 2 . The value of p will be chosen 
such that this correction will contribute to the order e 3//2 secular term and will give us the 
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required response. The equation of motion (|A1|) becomes 

he p 

x + uj 2 x = — ( e *M+ nT ) + c c ) _ e7i _ x 3 _ ^ A2 ^ 

and we try an expansion of the solution of the form 

x(t) = e 1/2 (A(T)e iu)t + c.c.) + e p+l ' 2 x p (t) + e 3/2 Xl (t) + ... (A3) 

Substituting this expansion into the equation of motion (|A2|) we obtain at order e 1 / 2 the 
linear equation as usual, and at order e p+1 / 2 

x p + oo 2 x p = I (Ae^ t+nT ^ + A*e inT + c.c.) . (A4) 

As expected, there is no secular term on the right-hand side, so we can immediately solve 
for x p , 

x p (t) = \ (-A_e^^) + ^ T + c.c) + 0(e). (A5) 

Substituting the solution for x p into the expansion ( |A3| ) , and the expansion back into the 
equation of motion contributes an additional term from the parametric driving which 

has the form 

e 2 P+ l/2^ ( _J_ e i{2u,t+nT) + ^tOT + \ r e i(ut+OT) + Q S 

4 \ 3w 2 u 2 J v ' 

h 2 (1 \ 
= e 2p+1/2 ^ ( 3^ + A*e i2nT J e iuJt + c.c. + non secular terms. (A6) 

To contribute to the order e 3 / 2 secular term, we see that we must set p = 1/2. This gives us 
the required contribution to the equation for the vanishing secular terms. All other terms 
remain as they were in Eq. (|T2|), so that the new equation for determining A(T) becomes 

2iU ^A _^_^ A + A *e l2nT ^j + icujA + 3\A\ 2 A + icur)\A\ 2 A = 0. (A7) 

Again, ignoring initial transients, and assuming that the nonlinear terms in the equation are 
sufficient to saturate the growth of the instability, we try a steady-state solution, this time 
of the form 

A{T) = ae inT . (A8) 
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FIG. 8: Response intensity | Ob I clS cl function of the frequency for fixed amplitude h = 1.5 
in the second instability tongue. Solid curves are stable solutions; dashed curves are unstable 
solutions. Thin curves show the response without non-linear damping (fj = 0). Thick curves show 
the response for finite nonlinear damping (fj = 1). 

The solution to the equation of motion (|7|) is therefore 

x (t) = e 1/2 (ae^ +en)t + c.c.) + 0(e), (A9) 

where the correction X\i2 of order e is given in Eq. ( |A"5| ) and, as before, we are not interested 
in the correction X\(t) of order e 3//2 , but rather in the fixed amplitude a of the lowest order 
term. We substitute the steady-state solution (|A8|) into the equation ( |A7| ) of the secular 
terms and obtain 



3\a\ 2 -2uQ- 



2 /i 2 



3 4w 2 



+ iu (7 + r]\a\ 2 ) 



Au 2 



(A10) 



We divide both sides of the last equation by •yu and define the rescaled variables: a 



a/^/TuJ, Cl = fj = LOT], and h = /i/2a/7^ 3 , in terms of which we obtain after taking the 
magnitude squared of both sides, in addition to the trivial solution a = 0, the non-trivial 
response 



2., 



3\a\ 2 - 2tt - -h 2 ) +(l+fj\a\ 2 ) 



2\ 2 



h 4 



(All) 



Fig. H shows the response intensity | (X I clS cl function of the frequency for fixed am- 
plitude h = 1.5 in the second instability tongue. The solution looks very similar to the 
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response shown in Fig. [I] for the first instability tongue, though we should point out two 
important differences. The first is that the orientation of the ellipse, indicated by the slope 
of the curves for fj — 0, is different. The slope here is 2/3 whereas for the first instability 
tongue the slope is 1/3. The second is the change in the definition of h. The lowest ampli- 
tude required for having an instability, is again on resonance (Q = 0), and its value is again 



h = 1, but now this implies that h = 2a/t^, or that h it scales as ^fy. This is consistent 



with the well known result (see, for example, Ref. [10]) that the minimal amplitude for the 
instability of the n th tongue scales as 7 1//n . 

APPENDIX B: EXPLICIT EQUATIONS FOR TWO AND THREE COUPLED 
OSCILLATORS 

1. Two coupled nonlinear oscillators 

For two coupled oscillators (N = 2) we have 

71 2ix 
Qi = 3> Q2 = y , (Bl) 

u\ = \- U\ u 2 2 = l- ^A 2 , (B2) 

and we choose the reference frequency u to be u 2 , so that S 2 — Q 2 = 0, r 2 = 1, S± — 
2(ui - uj 2 )he = 5 > 0, and n = u;i/cj 2 = r. For A < 1, S ~ A 2 /e7 and r ~ 1 + A 2 /2. The 
first mode is the symmetric one with xi(t) = x 2 (t) and the second mode is antisymmetric 
with Xi(t) = —x 2 (t). The equations (|32]) for the rescaled complex amplitudes a± and a 2 are 



(fi — S)ra 1 — i—a-i + — a i — j (l a i| 2 Oi + 2|a2| 2 «i + a 2 a*) 
3 

— —iff \r\ai\ 2 ai + 2r\a 2 \ 2 a\ + (2 — r)*! 2 ,^] = 0, (B3a) 

Qa 2 — i^a 2 + -J1CL2 — ^ (|a 2 | 2 a, 2 + 2|a! | 2 a 2 + a\a^) 
3 

— if} [9\a 2 \ 2 a 2 + 2\a 1 \ 2 a 2 + (2r - l)a\a* 2 ] = 0. (B3b) 



The two single-mode solution branches, having the general form of Eq. ( P5| ) and labeled 
Si and ^2 in Fig. f§, are easily obtained by setting a 2 or ai to zero in the coupled equations 
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above, respectively. This yields the analytical forms of these solutions which are 
St : [\\ai\ 2 - 2r(Q - 5) Y + r 2 (l + ^fY = h 2 



Si'- (^|a 2 | 2 -^] 2 +[l + ^|a 2 | 2X2 



h 2 . 



(B4a) 
(B4b) 



2. Three coupled nonlinear oscillators 



For three coupled oscillators (N = 3) we have 



7T 

qi = 4' 



7T 

92 =2' 



3tt 

g3 = T 



cu 2 = 1 - A 2 + A 2 /V2, u 2 2 = l- A 2 , u\ = 1 - A 2 - A 2 /V2, 



(B5) 
(B6) 



and we choose the reference frequency u to be a; 2 , so that <5 2 = 0, r 2 = 1, and 5% = 
2(ui - uj 2 )ht = -6 3 = 5 > 0. For A < 1, 5 ~ A 2 /V2e-f and r li3 ~ 1 ± A 2 ^^. The 
equations (p2|) for the rescaled complex amplitudes a 1; a 2 , and a 3 are 



(S2 - ajriOi - z riai H ft,a x 

3 

— (3|ai| 2 a! — |a 3 | 2 a 3 + 4|a 2 | 2 a 3 + 2a 2 a 3 + 4|a 2 | 2 ai + 2a\a\ + 6|a 3 | 2 ai + 3a 2 a^ — 2|ai| 2 a 3 — a 2 a 3 ) 



4 
—if] 

+ 



3-2y/2 o 2 _ 
3ri|ai| <i\ 

y/2-1 



V2 + 1 2 _ 2-V2 
: r 3 |a 3 | a 3 H 



1 



(2r 1 \a 2 \ 2 a 1 + (2 - ri)a 2 a*) 



i (2r 3 |ai| 2 a 3 + (2n - r 3 )a 2 a;) + - (2r 1 |a 3 | 2 a 1 + (2r 3 - r^a 2 ^) 



(B7a) 



f2a 2 — ia 2 + ha 2 — — (2|a 2 | 2 a 2 + 2|ai| 2 a 2 + a 2 a 2 + 2|a 3 | 2 a 2 + a 2 a 2 + a^a 2 a 3 + a\a* 2 a^ + aia 2 a 3 ) 



2 \/2 

|«2| 2 a2 H 5 — (2|ai| 2 a 2 + (2r a - l)a 2 a 2 ) 



+ ^?(2|a 3 | 2 a 2 + (2r 3 -l)a 2 a;) 

zs 2 + V2 2 + V2 T » 
(iZ + d)r 3 a 3 - z — - — r 3 a 3 H — ha 3 



(B7b) 



(3|a 3 | 2 a 3 — |ai| 2 ai + 4|a 2 | 2 a 3 + 2a 2 a 3 + 4|a 2 | 2 ai + 2a 2 a^ + 6|ai| 2 a 3 + 3a 2 a 3 — 2|a 3 | 2 ai — a 2 ai) 



—ir\ 



3 + 2^2 
4~ 

\/2 + l 



> I 12- , V^-l ,_ |2 _ , 2 + 

3r 3 |a 3 | a 3 H n|ai| a x H — 



(2r 3 |a 2 | 2 a 3 + (2-r 3 )a 2 a*) 



1 



(2ri|a 3 | 2 ai + (2r 3 - ri)a 3 a*) + - (2r 3 |ai| 2 a 3 + (2r 1 - r 3 )a 2 a 3 ) 



0. 



(B7c) 
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Only one single-mode solution of the form of Eq. (|33"D exists in the case of three oscillators, 
and involves the second mode. It is obtained by setting a\ = = in the coupled equations 
above. The analytical expression for this solution is 

S 2 : (3|a 2 | 2 -n) 2 +(l + r/|a 2 | 2 ) 2 = /i 2 . (B8) 
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